In [1]:
set.seed(999)
options(scipen = 9)
options(warn = -1) 
source("./environment/libraries.R")
knitr::opts_chunk$set(fig.height = 12, fig.width = 18, fig.dpi = 300)
knitr::opts_chunk$set(warning = FALSE)
All packages loaded successfully
In [2]:
name <- "Kenya_E1"
dataset <- read.csv(paste0("./test/", name, "_processed.csv"))
dataset_c <- data.frame(dataset[7:ncol(dataset)], row.names = dataset$Serial)
proj_coord <- data.frame("Easting" = dataset$Easting, "Northing" = dataset$Northing)
xy_coord <- data.frame("Longitude" = dataset$Longitude, "Latitude" = dataset$Latitude)

dataset_c_closed <- cbind(dataset_c, "Res" = 100 - rowSums(dataset_c)) 
head(dataset_c_closed)

dataset_spdf <- cbind(proj_coord, dataset_c_closed)
rownames(dataset_spdf) <- rownames(dataset_c)

structures_geojson_path <- file.path("./data/", paste0(name, "_topo_lines.geojson"))
structures <- st_read(structures_geojson_path, quiet = TRUE)
structures_utm <- st_transform(structures, crs = 32735) # Transform to UTM Zone 35S (WGS84)
A data.frame: 6 × 17
PCaTiMnFeNiZnRbSrYZrBaMgAlSiKRes
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
10.27313.47950.33090.06032.72480.00210.00330.00350.08080.00170.0123000000.09331.2275275.45467821.960071.97141662.32071
20.38353.95870.28280.05662.40540.00170.00310.00370.08390.00160.0041000000.09551.2341004.31907318.230192.22713966.70890
30.41584.07200.22650.06172.16340.00110.00330.00340.09210.00120.0045476330.07551.3168314.05277517.194472.48983467.82554
40.37043.95890.29020.06012.63740.00220.00350.00300.08370.00210.0102000000.08801.3842474.66847818.721382.42592465.29028
50.29813.39790.29840.07072.70030.00260.00300.00350.07970.00170.0152000000.10051.2800484.96270820.899372.48651463.39976
60.43734.00800.18080.05371.92700.00100.00290.00300.08600.00120.0035276360.05831.3045303.23393613.664892.97314372.06078
In [3]:
source("./utils/functions/create_quick_map.R")
options(repr.plot.width = 10, repr.plot.height = 10)

dbscan_result <- dbscan(dataset_spdf[, c("Easting", "Northing")], 
                        eps = 10, # Maximum distance between two points to be considered neighbors,
                        minPts = 2) # Minimum number of points required to form a dense region

Area <- factor(dbscan_result$cluster)
dataset$Area <- Area 
create_quick_map(dataset, structures, group_data = "Area")
No description has been provided for this image
In [4]:
dataset_spdf$Area <- Area 

# Create SpatialPointsDataFrames (compositional, ilr, clr)
spdf_comp <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                    data = dataset_spdf[, -c(1, 2, ncol(dataset_spdf))], 
                                    proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))
spdf_ilr <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                   data = as.data.frame(ilr(dataset_spdf[, -c(1, 2, ncol(dataset_spdf))])), 
                                   proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84")) 
spdf_clr <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                   data = as.data.frame(clr(dataset_spdf[, -c(1, 2, ncol(dataset_spdf))])),  
                                   proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84")) 
In [5]:
par(bg = "white")
options(repr.plot.width = 20, repr.plot.height = 15)
pairsmap(data = spdf_clr@data, loc = spdf_clr@coords)
No description has been provided for this image
In [6]:
par(bg = "white")
source("./utils/functions/create_grid_from_spdf.R")
source("./utils/functions/idw_compositional.R")
grid <- create_grid_from_spdf(spdf_ilr, resolution = 0.2, buffer = 2, convex_hull = FALSE)

# Perform IDW interpolation with compositional data
idw_compositional(spdf_ilr, orig = dataset_c_closed, grid = grid,
                  idp = 2, maxdist = 5, focal_window = 1,
                  plot_idw = TRUE, plot_ncol = 3)
No description has been provided for this image
In [7]:
par(bg = "white")
options(repr.plot.width = 20, repr.plot.height = 15)
swath(acomp(spdf_comp@data), spdf_comp@coords[,"Easting"], col = 8, xlab = "Easting", commonScale = FALSE)
swath(acomp(spdf_comp@data), spdf_comp@coords[,"Northing"], col = 8, xlab = "Northing", commonScale = FALSE)
No description has been provided for this image
No description has been provided for this image
In [8]:
par(bg = "white", mfrow = c(2, 1))
options(repr.plot.width = 15, repr.plot.height = 12)
source("./utils/functions/lag_distance_from_spdf.R")
source("./utils/functions/site_diagonal_from_spdf.R")
source("./utils/functions/create_gstat_from_spdf.R")


g <- create_gstat_from_spdf(spdf_clr, method = "universal") # Compute omnidirectional variograms
lag_dist <- lag_distance_from_spdf(spdf_clr) # Calculate the lag distance for variogram width
site_diag <- site_diagonal_from_spdf(spdf_clr) # Calculate site diagonal for variogram cutoff

v_omni <- variogram(g, 
               width = lag_dist / 2, cutoff = site_diag / 3, 
               cross = FALSE) 
plot(v_omni, type = "p")
No description has been provided for this image
No description has been provided for this image
In [9]:
par(bg = "white")
options(repr.plot.width = 18, repr.plot.height = 12)
source("./utils/functions/quick_anis_variogram.R")


g <- create_gstat_from_spdf(spdf_clr, method = "universal") # Create gstat object using CLR-transformed data

# Compute anisotropic variograms (alpha indicates the direction: 0 is NS, 45 is NE-SW, 90 is EW, 135 is NW-SE)
v_dir <- variogram(g, 
                   width = lag_dist / 2, cutoff = site_diag / 3,
                   alpha = c(0, 45, 90, 135), 
                   tol.hor = 22.5, 
                   cross = FALSE) 

quick_anis_variogram(v_dir, directions = c(0, 90), direction_labels = c("NS", "EW"))
quick_anis_variogram(v_dir, directions = c(45, 135), direction_labels = c("NE-SW", "NW-SE"))
No description has been provided for this image
No description has been provided for this image
In [10]:
par(bg = "white")
options(repr.plot.width = 15, repr.plot.height = 15)

# Create anisotropy plot with minimal parameters
vario <- logratioVariogram(data = acomp(spdf_comp@data), 
                          loc = spdf_comp@coords, 
                          maxdist = site_diag / 6,  
                          azimuth = c(0, 45, 90, 135), 
                          azimuth.tol = 45)  
image(vario)
No description has been provided for this image